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SUMMARY 


The  open-ended  coaxial  line  with  a  ground  plane  or  effectively  infinite 
flange  has  attracted  much  interest  in  recent  years  as  a  device  for 
non-destructive  measurement  of  complex  permittivity,  eg  in  bio-medical 
applications.  Existing  theoretical  treatments  of  this  configuration  all  involve 
approximations;  this  paper  presents  a  rigorous  treatment  using  a  variational 
approach. 
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THE  OPEN-ENDED  COAXIAL  LINE:  A  RIGOROUS  VARIATIONAL  TREATMENT 
T  E  Hodgetts 


INTRODUCTION 


This  paper  presents  a  rigorous  theoretical  analysis  of  the  open-ended  coaxial  line 
with  an  infinite  ground  plane  or  flange  4E*g",,TJ>_radiating  into  a  semi-infinite  lossy 
medium.  This  configuration  has  become  popular  for  measuring  the  complex  permittivity  of 
liquids  and  bio-medical  specimens,  because  it  does  not  require  a  sample  of  a  special 
shape.  However,  the  theoretical  analysis  has  previously  only  been  performed  by 

approximate  methods;  the  best  of  these  seems  to  be  the  point-matching  technique -j 
employed  by  Mosig  et  al  [1],  (discussed  with  attention  to  numerical  detail  bv^Syrrmr  (ij), 
but  there  are  several  others  in  the  literature  (see  the  bibliography  oj^  [3]).  GThe  present 
treatment  employs  the  variational  methods  described  by  Collin  [4Jand  /applied  to  the 
behaviour  of  enclosed  coaxial-line  discontinuities  by  Whinnery  et  al  [5],  Sornlo  [6],  Bianco 
et  al  [7]  and  the  present  author  [8].  a'-— i  /  .  ,  / 


THE  ELECTROMAGNETIC  FIELD  EQUATIONS  </ 


ha  <  Lcl^ 


We  begin  by  outlining  the  forms  of  the  fields  relevant  to  the  problem.  These  have 
been  discussed  at  length  by  many  writers  (e.g  Jones  [9],  Kerns  and  Beatty  [10]  and  the 
present  author  [8]),  so  this  section  will  be  little  more  than  a  statement  of  notation  without 
further  formal  references. 


A  region  of  space  which  is  linear,  isotropic  and  homogeneous,  and  contains  no  free 
charges  or  current-carrying  elements,  will  support  a  monochromatic  electromagnetic  field 
specified  by  the  six  space  components  of  the  conventional  complex  field  vectors  E  and  H, 
which  are  related  by  the  complex  monochromatic  forms  of  Maxwell's  equations: 


curl  E  +  juipH  -  0 


div  H  -  0 


curl  H  -  jo)e£  -  0 


div  E  -  0 


(1) 


In  eqns  (1),  o>  is  the  angular  frequency,  so  that  the  time  variation  of  the  fields  is  as 

exp(jut)  where  j  is  the  square  root  of  (-1),  and  p  and  r  are  the  permeability  and 
complex  permittivity  of  the  medium  in  the  space  region.  The  complex  permittivity  is 

defined  in  the  usual  way,  having  a  real  part  equal  to  the  ordinary  absolute  permittivity 
and  an  imaginary  part  equal  to  (— tr/cu)  where  a  is  the  conductivity  of  the  medium.  In 
our  problem  there  are  two  regions,  bounded  by  an  interface  between  them  (at  the 
end-plane  of  the  coaxial  line)  and  by  various  perfect  conductors  (the  flange  and  the 
cylinders  of  the  coaxial  line);  eqns  (1)  apply  separately  to  each  region. 

Both  regions  have  rotational  symmetry  about  the  axis  of  the  coaxial  line,  so  it  is 

convenient  to  use  cylindrical  polar  co-ordinates  (p,  <p,  z).  When  this  is  done,  it  turns 

out  that  the  six  space  components  of  E  and  H  de-couple  into  two  sets,  (H^,  E  _,  Ez) 
and  (E^,,Hp,Hz).  The  field  in  the  coaxial  line  far  from  its  end  necessarily  has 
components  belonging  only  to  the  first  of  these  sets  (at  normal  working  frequencies);  and 
it  is  possible  to  satisfy  all  the  conditions  at  the  interface  and  the  bounding  conductors 
without  introducing  components  in  the  second  set,  which  can  therefore  be  discarded. 

Equations  (1)  comprise  eight  scalar  equations.  Remembering  that  we  have  rotational 
symmetry  -  so  that  all  the  field  components  are  independent  of  <p  -  and  discarding  the 
four  equations  which  relate  to  the  irrelevant  (null)  components  (Ev„,  Hp,  Hz),  we  have 
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four  significant  equations,  one  of  which  is  an  immediate  consequence  of  the  other  three. 
(This  redundancy  follows  directly  from  eqns  (1);  if  we  take  the  divergences  of  the  two 
(vector)  equations  involving  curls,  the  other  two  (scalar)  equations  follow  immediately.) 
The  three  remaining  equations  are 


E 

z 


(PV  -  \  h 


(PH  )  +  - y 

*  8? 


(PH  )  +  k  (pH  ) 
< fi  V 


-  0 


(2) 


(3) 


where  k2  =  (as  usual)  and  the  unit  vectors  (£,  g,  z)  form  a  right-handed  set  (see 

Fig  1). 

MODES  IN  THE  COAXIAL  LINE 

The  tangential  electric  field  on  any  perfectly-conducting  surface  is  zero.  The 
directions  tangential  to  the  cylinders  of  the  coaxial  line  are  g  and  z;  hence,  from  eqns 
(2), 


8 

3? 


(pH  ) 


0  at  p  -  r  and  p  -  R 


(4) 


where  r  and  R  are  the  radii  of  the  inner  and  outer  conductors. 

This  condition  has  no  explicit  z-dependence,  so  eqn  (3)  can  be  solved  subject  to  it 
by  separating  the  variables  p  and  z.  The  equation  in  z  is  harmonic  or  exponential;  the 
equation  in  p  is  a  close  relative  of  Bessel's  equation.  A  separation  constant  appears, 
which  is  constrained  to  take  only  particular  values  (by  eqn  (4)).  The  complete  field  for 
z  <  0  (within  the  line;  see  Fig  1)  can  now  be  written  down  (after  [8]  and  [9])  as 
follows. 

The  equation 


J0(kir)Y0(kiR)  ‘  VkiR)Y0(kir>  “  0 


(5) 


has  a  set  of  real,  positive  and  simple  zeros  kj,  where  i  is  an  integer  numbering  the  zeros 
in  ascending  magnitude,  and  Jq  and  Yg  are  Bessel's  functions  of  the  first  and  second 
kinds  of  order  0  (Watson  [11]).  For  each  i  we  define  a  constant  Cj  satisfying 


Cj  -  J^kjrj/Yofkjr)  -  J0(kjR)/Y0(kjR) 


(6) 


(which  is  consistent  with  eqn  (5)),  and  we  also  define  three  functions  of  p,  thus 
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W) 


W>  +  CiYQ(kip) 


(7) 


where  a  is  one  of  the  integers  0,  1  or  2  (2  being  only  rarely  required);  the  ZQ  will  only 
be  used  with  arguments  of  the  form  (kjp)  so  that  the  required  value  of  Cj  is  always 
obvious.  From  eqns  (6)  and  (7),  the  mode  constants  kj  now  satisfy 


Z0(kir)  -  ZQ ( k  j  R )  -  0 


(8) 


for  all  i.  Also,  from  [8]  or  [9], 


^  (xJ  |  (x) )  -  xJq(x)  and  ^  (xYj(x))  -  xYQ(x), 


so 


d(k.p) 


[(k.p)Z1(k1p)]  -  (kiP)Z0(k.P) 


(9) 


The  field  components  then  are 


H  (z  <  0)  -  —  art 

»?AP  o 


-jRAZ  r  JkA2 

!  -  PAe 


♦  y 

i-1 


yiz 

a.e  ZjCk.p) 


(10a) 


E  (z  <  0) 
P 


1 

p  ao 


r  -JkAz  r  jkAz 

+  rA  e 


V  7iz 

+  )  a^  ZjCkjP) 


1-1 


(10b) 


Ez(z  <  0) 


.  \ 

L 

1-1 


ki  1 


ai  e  Z0(kip) 


(10c) 


where  the  suffix  A  denotes  region  A  in  Fig  1,  and 


2  2  2 
kA  ’  "  VA  and  7i 


k?  -  k»  and  rj. 

1  A  A 


(ID 


with  k^,  tja  and  the  7;  being  positive  (if  is  purely  real),  or  having  positive  real  parts 
(if  is  complex).  (It  is  assumed  that  the  frequency  is  low  enough  for  all  the  yf  to 
have  positive  real  parts,  so  that  only  the  ordinary  coaxial  transmission-line  mode 
propagates  freely,  all  other  modes  decaying  exponentially  away  from  the  plane  z  =  0  even 
in  the  lossless  case.)  The  remaining  parameters,  aQ  and  the  aj,  are  to  be  determined 
later.  This  representation  may  be  compared  with  that  in  [2]. 

THE  FIELD  IN  THE  UNBOUNDED  DIELECTRIC  REGION 

The  specification  of  the  unbounded  field  is  considerably  more  complicated.  Symm 
[2]  takes  it  directly  from  a  result  given  by  Lewin  [12];  however,  as  we  have  started  from 
first  principles  (eqns  (1)),  it  is  appropriate  to  take  the  treatment  somewhat  further  back  so 
as  to  show  more  clearly  the  connection  between  the  results  for  the  bounded  and 
unbounded  cases,  and  to  justify  the  transformation  of  variables  used  by  Symm  [2]. 

We  begin  with  a  theorem  quoted  in  [12)  from  Baker  and  Copson  [13];  the  proof  of 
this  is  an  exercise  in  pure  mathematics,  not  physics,  and  the  reader  is  referred  to  [13]  for 
the  details.  (To  apply  the  theorem  we  change  temporarily  to  Cartesian  co-ordinates, 
retaining  the  z-direction  and  taking  the  x-z  plane  as  the  p  =  0  plane;  see  Fig  1). 

Let  i'  be  a  solution  of  Helmholtz'  equation  in  region  B  of  Fig  1;  that  is,  V^|.  + 
kg^,,  s  o,  where  kg^  =  oj^gfg  and  kg  is  positive  (for  real  tg)  or  has  a  positive  real 
part  and  a  negative  imaginary  one  (for  complex  eg)  (compare  eqns  (11)  and  the  text 
immediately  following  them).  Let  v  further  have  no  singularities  in  the  positive-z 
half-space,  and  satisfy  the  radiation  condition  at  infinity;  and  let  the  function  (8r/dz)(z=Q) 
(defined,  if  necessary,  as  the  limit  as  z  tends  to  zero  from  above)  be  denoted  by  t!(x,y). 
Then,  in  the  positive-z  half-space, 


.■(x'  ,  y'  ,  z’ ) 


exp(-jk  D) 

J  J  ^(x,y)  - p -  dxdy 


(12) 


where  D  is  the  positive  square  root  of  ((x  -  x')2  +  (y  -  y')2  +  z*2)  anc|  the  double 
integral  extends  over  the  whole  z  =  0  plane.  (As  before,  a  harmonic  time  factor  exp(jwt) 
has  been  suppressed.) 

Now,  region  B  satisfies  the  same  symmetry  conditions  as  region  A,  from  which  it 
was  deduced  that  the  field  in  region  A  has  only  components  ( H Ep,  Ez).  We  infer 
that  the  same  is  true  of  the  field  in  region  B,  since  the  two  regions  are  connected  by 
continuity  conditions  on  the  tangential  components  of  E  and  H  at  the  z  =  0  plane.  We 
cannot,  however,  apply  eqn  (12)  immediately  to  the  reduced  set  of  field  components, 
because  the  critical  components  H^  and  E„  do  not  satisfy  Helmholtz*  equation;  but  the 
Cartesian  components  of  E  and  H  do  satisfy  it.  It  is  easiest  to  establish  this  result  and 
deduce  its  consequences  if  we  follow  the  treatment  in  [9],  where  it  is  shown  that,  subject 
to  the  assumptions  already  made,  the  electromagnetic  field  in  a  region  such  as  B  can  be 
expressed  in  terms  of  two  scalar  functions  satisfying  the  damped  wave  equation  (normally 
taken  as  the  magnitudes  of  Hertzian  vectors  with  directions  fixed  along  the  z-axis).  Using 
the  implied  time-harmonic  factor  with  the  results  in  [9],  the  Cartesian  field  components 
are 


_  a  n  .  8m  u  3  m  .  an 

Ex  “  ^  *bhx  ”  353Z  +  jwVb 
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(13) 


a2n  ^  .  3m 
3^35  +  JW  dt 


^BHy 


92m  an 

3^35  '  Ja,VB  3x 


E 


^  +  R2n 

az2  B 


"BHz 


32m  ^  .2 

a?  B 


where  the  complex  scalars  n  and  M  now  satisfy  Helmholtz'  equation. 
We  now  wish  to  prove  that 


dH 

y 

3 T 


dH 


-  J“fBEx  and  3T  “  jufBEv 


(14) 


Given  that  Hz  is  zero  in  region  B,  these  results  in  fact  follow  directly  from  eqns  (13). 
A  less  direct,  but  more  physically  instructive,  way  of  establishing  eqns  (14)  is  by 
transforming  into  the  (p,  p,  z)  co-ordinates  and  discarding  derivatives  of  n  and  M  with 
respect  to  p  -  the  E  and  Ji  fields  are  independent  of  p,  so  common  sense  would  be 
contradicted  if  n  and  M  were  not  also  independent  of  it.  We  then  find  that  the 
components  (H^,  Ep,  Ez)  depend  only  on  FI,  and  (E^,,  Hp,  Hz)  depend  only  on  M;  as 
this  second  set  of  components  is  null,  we  may  discard  M,  and  eqns  (13)  then  simplify  to 
forms  from  which  eqns  (14)  follow  without  further  assumptions. 

It  is  also  clear  from  eqns  (13),  whether  or  not  the  terms  in  M  are  discarded,  that 
Hx  and  Hy  satisfy  Helmholtz'  equation  -  because  FT  and  M  do  so,  and  therefore  so  do  all 
their  Cartesian  derivatives  (as  can  be  shown  by  differentiating  Helmholtz'  equation  and 
inverting  the  order).  Hence  we  may  set  Hx  or  Hy  for  v  in  eqn  (12),  taking  the 
corresponding  (dv/dz)  from  eqns  (14);  then 


H  (x' , y'  , z ' ) 


e*P(-jkBD) 

Ex(x,y,0)  - g -  dxdy 


Hx(x' ,y' ,z') 


exp(-jkBD) 

E  (x,  y ,  0)  - g -  dxdy 


(15) 


Now,  by  taking  components  in  the  various  orthogonal  directions  in  Fig  1  it  is  easy  to 
show  that 


H  -  H  cos  if>'  -  H  sin  <p'  (in  the  primed  co-ordinates) 

v  y  x 


and 


E 

y 


E  cos  -  E  sin  ip 
P  v° 


E  s i n  p  +  E  cos  p 
p  p 


(in  the  unprimed  co-ordinates) 
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using  x  =  p  cos  <p,  y  =  p  sin  ip  and  similar  primed  equations.  We  substitute  from  these 
equations  for  Ex  and  Ey  in  eqns  (15);  then  we  multiply  the  first  of  eqns  (15)  by  (cos  vV) 
and  the  second  by  (-  sin  ip'),  add,  and  combine  the  two  integrals  into  one.  We  thus 
obtain,  remembering  that  E^  is  zero, 


H  (x*  ,y’ ,z’ ) 


exp(-jkBD) 

Ep(x,y,0)  cos  (p  -  <f>')  - g -  dxdy 


We  now  complete  the  conversion  back  to  cylindrical  polar  co-ordinates.  Since  Ep  is 
tangential  to  the  z  =  0  plane  over  which  we  are  integrating,  it  vanishes  for  all  radii  p 
less  than  r  or  greater  than  R,  since  these  parts  of  z  =  0  are  occupied  by  perfect 
conductors  (Fig  1);  and  the  differential  area  (dxdy)  transforms  into  (pdpdp),  as  is  well 
known.  The  last  equation  then  becomes 


H  (P’  ,p’  ,Z’  ) 
<f> 


R  2ir 

exp(- jkgD) 

E^(p,(p,0)cos(lp  -  p' )  - g -  pdpdp  (16) 

r  0 


where  now  D 2  =  p2  +  p’2  _  2pp*  cos(<p  -  p')  +  z ’2,  and  the  dependence  of  Ep  on  <p 
and  of  on  <p'  is  only  formal.  This  is  the  result  given  by  Symm  [2], 

CURRENT,  VOLTAGE,  ADMITTANCE  AND  REFLECTION  COEFFICIENT 

It  is  possible  to  carry  on  the  analysis  from  here  in  terms  of  the  principal  mode's 
"voltage  reflection  coefficient"  (introducted  in  eqns  (10))  or  in  terms  of  current, 
voltage  and  admittance.  The  latter  representation  is  quite  logical  (since  we  have  a 
terminated  two-conductor  system)  and  is  more  satisfying  to  microwave  engineers;  but  it 
requires  a  lengthy  discussion  to  justify  it  rigorously  (eg.  [10]).  For  convenience  we  shall 
proceed  in  terms  of  the  reflection  coefficient;  the  link  with  admittance  can  then  be  made 
using  the  standard  equation  ([8], [10]) 


Y 

Y 


1  -  r. 


i  +  r. 


(17) 


where  Y0  is  the  characteristic  admittance  of  the  coaxial  line  and  is  2W(i}A2n(R/r)). 

INTEGRAL  RELATIONS  AND  ORTHOGONALITY 

It  is  well  known  that  the  modes  in  modal  expansions  like  eqns.  (10)  satisfy  integral 
relations  which  make  them  mutually  orthogonal.  The  necessary  relations  for  our  problem 
are  quoted  in  [8]  from  [11]  and  are 


}  [  -  )  Pdp  -  Cnp 

J  Zj  (0p)pdp  -  \  p2[Zj(flp)  -  Zo(0p)Z2((?p)] 
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f  [lw>p)pip  -  -  Jv">> 


I  Zj (»1P)Z1 («2p)pdp 


[fl2pZ1(^p)Z0(e2p>  -  «1pz1(^p>z0(e1p)] 


(18) 


(0J  *  «2) 


<*1  - 


The  last  two  of  these  are  the  orthogonality  relations  properly  so  called,  as  can  be  seen  by 
taking  the  integration  from  r  to  R,  setting  the  0's  equal  to  one  or  two  of  the  kj  and 
using  eqn  (8). 

VARIATIONAL  FIELD  REPRESENTATIONS 

We  are  now  ready  to  tackle  the  heart  of  the  problem.  Let  the  field  components 
and  Ep  at  z  =  0  be  denoted  by  the  forms  H ^  and  Ep.  By  allowing  z  to  tend  to  0 
in  eqns  (10)  and  (16)  we  obtain  two  representations  of  H^\  these  must  be  equal,  since 
the  tangential  components  of  H  are  continuous  across  the  plane  z  =  0,  but  in  some  of 
what  follows  it  will  be  convenient  to  distinguish  them,  as  H+  and  H~  according  to  the 
sign  of  z.  and  Ep  are  functions  only  of  p,  because  of  the  cylindrical  symmetry. 

From  eqns  (10)  and  (16),  with  some  interchanges  of  p  and  p'  for  convenience, 


•ft*1  -  rA>  “  VAP 


*  1-1 


J««i 


aiZl (kjP) 


V 


*  i-1 


i  J 


a.Zj (k.p) 


V 


V 

i-1 


J“«, 


J 


ai^l (kjP) 


(19) 


R  2tt 


Ep(p' ) cos (<p  -  <p' )  D 


exP(~ jkgDQ) 


p'd^'dp' 


r  0 


where  D„  - 


2  '2 

p  +  p  -  2pp'cos(ip  -  ip' ) 


It  is  possible  for  the  integrand  in  eqn  (19)  to  become  singular;  to  deal  with  this  problem 
a  suitable  approach  to  the  limit  will  be  used. 

Now,  from  eqns  (10),  (18)  and  (8),  with  interchanges  of  integration  and  summation. 
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aQ(l  +  rA)Cn(R/r) 


K  oo 

“  1  K  '  I  ai2i(kip)]  dp 


i-i 


£  R 


’  I  Ep  dp  "  l  ai  I  Zl(kip)dp 

r  i-1  r 

R 

-  J  E  dp  ;  and  similarly 


J  Z2(kiP)pdp 


-  |  Zj (kjp)p 


a0d+rA)  S 

£p-— jr^-l  WVP) 

i'-i 
i  '*i 


dp 


-  j  EpZ1(kip)pdp  . 
r 

Substituting  these  results  for  ag  and  aj  into  eqn  (19)  gives  us 


2r 


r?ACn(R/r) 


i  -  r. 


i~rr 


I  V’ 


2rp 


00 

V 

R 

f  2 

-1 

R 

> 

i-l 

*i  J 

Zj(k.p) 

J  Z1(k.p)pdp 

r  J 

1 

r 

E  Z. (k  p)pdp 
P  1  > 


JWf, 


R  It 


I  \  Ep('P'  )C0S(-'r'’  ~  ) 

r  0 


exp(-jkBDQ) 


p'dvp'dp’ 


This  equation  contains  the  quantity  (2*7(T7ACn(R/r)))((l  -  T^/fl  +  P /^)),  which  by 
(17)  we  are  entitled  to  call  the  terminating  admittance  Y  of  the  coaxial  line.  Let  us 
multiply  the  equation  by  an  arbitrary  function  Tp  of  p  and  integrate  from  r  to  R; 
then  have 
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>fA 

R 

f  2 

-i 

R 

r 

R 

r  — 

>i 

j  Z^ (kj p)pdp 
r 

j  E^ZjfkjPlpdp 
r 

J  EpZj (kjP)pdp 
r 

eqn 

now 

we 
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*  R  2.t_  rexp(-jk_D0) 

+  J"*B  J  J  I  Ep(p)Ep(p,)cos(^  -  *'>  - 67 - 

r  r  0  0 


pp’  dip’  dp '  dp 


In  eqn  (20)  the  function  Ep  is  the  exact  form  of  the_  radial  electric  field  Ep  at  the 
plane  z  =  0.  We  desire  now  to  choose  the  function  Ep  so  that  we  obtain  the  best 
possible  value  for  Y  when  we  use  (as  we  inevitably  must)  an  approximation  to  Ep.  To 
do  this,  we  consider  the  effect  on  eqn  (20)  of  a  small  change  in  Ep  or  E p.  In  the 
terminology_of  the  calculus  of  variations  (eg  [4],  [8])  the  first  variation  of  eqn  (20)  with 
respect  to  Ep  is 


R  R 

I  Vp  I 

r  r 


&E  dp  +  6Y 
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R  R 

f  Vp  1 

r  r 


Y  r jci)<  i  R  i-l  r  r  1  r  R  _ 

*  2  — —  J  Z1(kjp)pdp  |  EpZ1(k.p)pdp  \6EpZi  (k.p)pdp 

i-l  *■  t  J  l.  r  JLr  -ILr 


R  R  It 


+  jw«B  in  5E  (p)£  (p')cos(p  -  p') 


r  r  0 


exp(-jkBDQ) 


pp ' dp' dp ' dp 


or  »y(  |  £dp][  |C  dp]  - 


R  _  f  “  rjuci  f  R  2  1-1  R 

|  dp6£  2ip  )  — — —  Zj (k. p)  J  Zj(k.p)pdp  J  £  Zj (k.p)pdp 

r  ,ri l  Ti  J  r  r 


+  JWfBP 


R.  2tt  rexp(-jk  D0)‘ 

E  (p')cos(p  -  p’)  - -  p'dp'dp' 

r  0  0 


R 

'  ( Vp 

r 


The  integrand  in  brace  brackets  here  vanishes,  on  account  of  the  equation  preceding  eqn 
(20);  hence  6Y  is  zero  for  variations  of  Fp.  From  the  symmetry  of  eqn  (20),  we  can 
immediately  deduce  a  counterpart  to  eqn  (21)  with  £„  interchanged  with  Fp  (and  p 
interchanged  with  p'),  and  it  then  follows  that  setting  Fp  =  Ep  makes  6Y  vanish  for 
variations  of  Ep  also.  This  is  the  optimum  form,  since  first-order  errors  in  Ep  then  only 
cause  second-order  errors  in  Y.  The  vanishing  variation  property  can  also  be  used  in 
reverse,  after  Ritz  ([4], [8]),  as  follows.  Putting  Ep  =  Ep  in  eqn  (20)  gives  us 


R 

I  Vc 

i" 
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If  an  approximate  form  for  Ep,  containing  adjustable  parameters,  is  inserted  into  eqn  (22), 
it  will  represent  Ep  as  well  as  it  can  when  the  first  variations  of  the  equation  with 
respect  to  the  adjustable  parameters  all  vanish.  Before  applying  this  technique,  however, 
we  will  first  transform  the  triple  integral  in  the  equation  into  a  more  tractable  form. 

REPRESENTATION  OF  THE  PROPAGATION  FACTOR 

The  factor  (exp(-jkgDg)/Dg)  -  sometimes  called  the  propagation  factor  -  usually 
appears  in  analyses  of  this  kind,  and  has  a  large  body  of  literature  associated  with  it  in 
classical  electromagnetics  and  optical  diffraction  theory,  including  the  Sommerfeld  dipole 
problem  which  is  perhaps  the  present  problem's  nearest  relative  (see  Banos  [14]  for  a 
very  full  discussion  and  bibliography).  The  integral  transform  we  need  is,  in  fact, 
discussed  by  Bateman  [15]  in  the  context  of  Sommerfeld ’s  work.  However,  it  is  unsafe 
simply  to  quote  the  result  from  [15],  as  a  branch-point  integral  is  involved  and 
considerable  care  is  needed.  Accordingly,  we  shall  prove  it  directly,  following  the  method 
given  by  Jeffreys  and  Jeffreys  [16]. 

We  note  first  that  (exp(-jkgDg)/Dg)  is  the  limit  as  z'  tends  to  zero  of 
(exp(-jkgD)/D)  (see  eqns  (12)  and  (16)),  and  it  is  this  latter  form  which  we  shall 
represent,  so  that  the  limit  can  be  taken  later  (see  the  comment  following  eqn  (19)). 

Now, 


4.  h 

»xp(-0z' )Jo[Doy(kg  +  e2)j  “  57  {  exp["  c(z'  +  JDqCosq)  +  jkfiD0sinQ  j da 


where  Dg,  z'  and  kg  have  the  meanings  already  used,  C  is  a  general  complex  number 
with  a  non-negative  real  part,  and  Jg  is  the  Bessel's  function  of  the  first  kind  and  order 
0  already  introduced.  (The  sign  of  the  square  root  on  the  left-hand  side  is  of  no 
consequence,  since  Jg  is  an  even  function.)  This  result  is  easily  proved  by  removing  the 
factor  exp(-Gz’)  from  both  sides,  expanding  the  remaining  exponential  integrand  as  an 
absolutely  convergent  series  with  infinite  convergence  radius,  and  integrating  it  term  by 
term;  the  equation  then  reduces  to 


which  is  true  by  definition  (eg.  [16]). 
Next,  we  consider  the  integral 
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00 


JQ(sD0)sds  ■  I  (say) 


(24) 


exp(-z' (s2  -  k2) ^) 


(s 


w 

V 


in  which  the  path  of  integration  is  the  positive  real  axis,  and  we  introduce  the  temporary 
requirements  that  z'  is  strictly  positive  and  kg  is  strictly  complex;  the  square  root  has  its 
natural  definition  (approximately  equal  to  s  for  large  s  on  the  integration  path).  If  we 
change  the  variable  of  integration  to  this  natural  square  root,  which  we  will  call  C,  we 
then  have 


1 


00 


-  I 

jkg 


exp(-  z'C) 

- S - 


.CdC 


where  the  lower  limit  of  integration  is  fixed  by  continuity  (since  (s^  -  kg2)  has  a  positive 
imaginary  part  for  all  real  s)  and  the  square  root  in  the  argument  of  Jq  causes  no  trouble 
as  remarked  above.  Since  also  the  6  in  the  denominator  is  cancelled  by  the  C  in  the 
numerator,  we  have  an  integrand  which  is  well-behaved  in  the  entire  complex  finite 
C-plane;  this  makes  it  unnecessary  to  specify  the  path  of  integration  in  detail,  and  (more 
importantly)  allows  us  to  substitute  from  eqn  (23)  and  interchange  the  order  of  integration 
of  the  resulting  double  integral,  which  also  has  a  well-behaved  integrand;  so 


I 


1 

27 


2t  CO 

j  \  exp (~C (z' 
0  jkB 


+  jDQcosa)  +  jkgDQsinQr 


The  integration  over  Q.  can  be  performed  at  sight,  remembering  that  z'  is  strictly  positive 
so  that  exp(-z'G)  tends  to  zero  as  Q  tends  to  infinity;  we  then  have 


I 


1 

27 


2r 


\ 


exp(- jkgZ'  +  kgDgCOSQ  +  jk^sina) 

z~'  +  jD^cosa 


da 


This  integral  is  transformed  by  the  substitution  0  =  exp(ja),  which  (contrary  to 

appearances)  does  not  introduce  a  branch  point  (because  all  the  occurrences  of  a  are 
trigonometric,  except  the  differential  da  which  can  be  written  as  d0/(j  13)).  We  then  have 


I 


1 


(2ir) 

I 

(0) 


exp(-jk  z'  +  k  D  (3) 

- - - ^ -  d0 

0z'  +  i jDQ(p  +  1) 


where  the  path  of  integration  is  anti-clockwise  around  the  unit  circle  in  the  complex 
0-plane.  This  integrand  has  two  simple  poles  at  the  points  (jz’/Do)  *  j/(l  +  (z'/Dq)^), 
and  since  z’  is  strictly  positive  it  is  clear  that  the  contour  only  encloses  the  pole  at  0  = 
(jz’/Do)  -  j  y(l  +  (z'/Dq)^).  We  then  obtain,  from  eqn  (24), 
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00 


Jq(sDq)s<1s 


I 


(25) 


J 


exp(-z' (s2  -  k2) *) 


(s' 


exP(” jkgD) 
D 


remembering  that  =  Dq^  +  z'2  from  the  definitions  of  D  and  Dp.  This  form  may 
now  be  used  in  the  triple  integral  in  eqn  (22),  with  the  limit  as  z'  tends  to  zero  being 
taken  later  -  this  is  physically  justified  by  the  continuity  of  the  transverse  field 
components  at  the  plane  z  =  0,  and  is  logically  justified  by  the  flow  of  the  argument 
from  eqn  (16)  to  eqn  (22).  The  triple  integral  accordingly  becomes  the  limit  as  z'  tends 
to  zero  of 


I 


r 


I 


r 


E (p)E (p')pp' 

r  r 


2^ 

|  cos(<p 
0 


“  exp(-z'(s2  -  k2)*) 

**>  f  -  ;-2  ,2^  VsDo)sds  d*' 

0  (s  -  V 


dp'  dp 


INTEGRATION  OVER  p' 

We  now  consider  the  term  in  square  brackets  in  the  last  expression,  still  maintaining 
the  temporary  requirements  that  z'  is  strictly  positive  and  kg  is  strictly  complex.  As 
before,  we  can  remove  the  branch-point  by  temporarily  transforming  from  the  variable  s 
to  the  previously  defined  variable  6.  We  may  then  invert  the  order  of  integration  and 
transform  back  to  the  variable  s,  and  the  bracketed  part  of  the  expression  becomes 


00 


I 


exp 


2t 

J  J0(sD0)cos(v? 

0 


)dp' 


ds 


(26) 


The  integral  in  brace  brackets  can  be  evaluated  in  closed  form;  and,  remarkably,  the 
integration  variables  p  and  p'  are  separated  in  the  process.  We  remark  first  that  Dp 
depends  on  < p'  only  through  cos(p  -  p'),  and  this  property  is  preserved  for  the  whole 
integrand;  and,  moreover,  if  F  represents  an  arbitrary  well-behaved  function  we  have 


It  It 

J  F^cos(p  -  p')jdp'  -  |  F^cos(p'  -  p)jdp' 

0  0 

(2  *-p) 

-  j  F(cos  p’ )  dy>' 

-s p 

(2ir-p)  0 

-  |  F(cos  p’  )  dp'  +  J  f|cos(2jt  +  p')jd(2r  +  p' ) 

0  -ip 
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(2  T~p)  2  IT 

-  |  F(cos  <p')d<p'  +  J  F(cos  if>' ) dy?' 

0  (2*-<p) 

2ir 

-  |  F(cos  )d<p' 

0 


Hence  the  brace-bracketed  integral  may  be  written  as 


2r  _ _ _ 

|  J0[/(sp)2  +  (sp')2  -  2(sp) (sp'  )cosi/>' j  costs' di^' 
0 


using  the  definition  of  Dq.  Now  Jq  has  a  triangular  expansion  given  in  [16]  and  [11]  as 
(in  our  notation) 


y^7sp7^-+-(sp^)^— ^-^(^pHsp^cosp' 


-  J0(sp)J0(sp' )  +y  2005(1^' )Jm(sp)Jm(sp' ) 
m-1 

where  Jm  denotes  the  Bessel  function  of  the  first  kind  and  of  order  m.  This  series  has 
an  exponential  tail,  since  the  asymptotic  formula  in  [9]  for  the  Bessel  functions  of  the 
first  kind  of  fixed  argument  and  large  variable  order  shows  that  Jm(w)  (regarded  as  a 
function  of  m)  decays  exponentially  to  zero  as  soon  as  m  significantly  exceeds  ($ew),  e 
being  the  base  of  natural  logarithms.  We  may  therefore  interchange  the  order  of 
integration  and  summation  after  substituting  this  series  into  the  expression  above  for  the 
brace-bracketed  integral  in  eqn  (26).  All  the  integrals  over  <p'  then  vanish,  except  the 
one  where  the  index  m  is  1 ;  and  the  brace-bracketed  integral  is  simply  given  by 


{  }  -  2»J1(sp)J1(spt) 

Hence,  from  eqn  (26)  and  the  expression  before  it,  the  triple  integral  in  eqn  (22)  is  the 
limit  as  z'  tends  to  zero  of 


exp(-z’(s2  -  k2)*) 

- x - *— J -  2irsJ.  ( sp )  J .  (sp' ) dsdp ' dp 

( ss  -  k‘)' 


R  R 

J  J 


E  (p)£  (p’)pp' 
P  P 
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(This  method  of  integrating  over  <p'  is  described  by  Morse  and  Feshbach  [17]  in  the 
context  of  a  similar  problem.) 

THE  RITZ  PROCEDURE 

We  now  apply  to  the  exact  equation  (22)  the  Ritz  procedure  of  successive 
approximations  described  immediately  after  that  equation.  We  want  an  approximate  form 
for  Ep  which  depends  in  a  convenient  way  on  adjustable  parameters,  approaches  the  exact 
form  if  enough  parameters  are  used,  and  (for  stability)  satisfies  the  boundary  and 
consistency  conditions  (implied  by  eqns  (2)  and  (4))  at  every  approximation  stage.  The 
natural  form  is  the  one  obtained  by  truncating  the  series  expression  for  Ep  at  z  =  0  from 
eqn  (10),  namely 


ECp)  -  i  «n  +  7  k  ®  Z  (k  p)  (27) 

p  p  U  L  m  m  l  m 

m— 1 

for  the  N1*1  approximation,  where  oq  and  the  am  are  the  adjustable  parameters  (the 
factor  km  has  been  included  for  later  convenience).  The  mode  functions  1/p  and  the 
Zi(kmp)  are  an  orthogonal  complete  set  (see  [4]  and  [8]),  so  the  form  (27)  can  represent 
as  accurately  as  we  please  in  the  mean  any  function  in  the  relevant  range  of  p  and 
satisfying  the  relevant  boundary  conditions,  if  N  is  sufficiently  large.  At  each 
approximation  stage  we  substitute  £  p  into  eqn  (22),  adjust  every  parameter  for  a  vanishing 
first  variation,  and  calculate  the  corresponding  value  of  Y  (say,  Yjq  for  the  Nth  order). 
Since  £p  depends  linearly  on  its  adjustable  parameters,  there  is  at  each  stage  only  one 
stationary  point  where  the  variations  with  respect  to  every  parameter  vanish;  hence,  since 
the  exact  £p  is  also  stationary  in  this  sense,  as  N  tends  to  infinity  £p  tends  to  £„  in  the 
nunn  and  Y^  tends  to  Y,  and  so  Y  can  be  found  by  calculating  suitable  Yjq  and 
extrapolating  to  infinity  using  a  suitable  procedure  for  accelerating  convergence. 

We  substitute  from  eqn  (27)  first  into  the  expression  above  for  the  triple  integral  in 
eqn  (22).  The  result  is  a  finite  double  sum  of  triple  integrals,  and  the  limit  as  z'  tends 
to  zero  may  be  taken  subsequently  in  each  of  these  (provided  they  all  remain  convergent, 
which  we  shall  see  they  do);  we  then  have  for  this  part  of  eqn  (22) 
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exp(-z' (s  -  kg) J) 
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where  the  integrals  I  satisfy  Imn  =  Inm  and  are  given  by 
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exp(-z'(s2  -  k2)*) 

- x - x — j -  2tsJ.  (sp)J1  (sp' )dsdp'dp  (28a) 
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I 

mn 
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I  I 


r  r 


“  exp(-z' (s2  -  k?) b 

I  -— r— O - (V>Z1(V><V>Z1(V>- 

0  (s  -  kg) 2 


.  2irsJj  (sp)Jj  (sp‘  )dsdp'dp  (28c) 


Still  maintaining  the  temporary  requirements  that  z'  is  strictly  positive  and  kg  is  strictly 
complex,  we  can  once  more  remove  the  branch-point  from  each  of  the  I-integrals  by 
temporarily  transforming  from  the  variable  s  to  the  previously-defined  variable  C;  we  may 
then  invert  the  order  of  integration,  transform  back  to  the  variable  s,  and  integrate  first 
over  p'  and  p.  The  p'  and  p  integrals  are  all  special  cases  of  eqns  (18),  with  a 
Z-function  which  is  a  degenerate  form  with  Cj  (in  the  notation  of  eqn  (7))  always  zero. 
As  a  consequence  of  the  fixed  Cj  it  does  not  matter  that  the  arguments  of  J}  are  not  of 
the  form  (kjp)  or  (k;p'),  and  on  applying  eqns  (18)  and  recalling  eqn  (8)  we  obtain 
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|  J  j (sp)dp  -  s_1[j0(sr)  -  Jq(sR)] 

r 

(29) 
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f  (k  p)Z1 (k  p)J. (sp)dp 
J  m  i  m  l 
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(k  r)Z.(k  r)sJ.(sr)  -  (k  R)Z.(k  R)sJn(sR) 
m  i  m  u  m  i  m  u 


with  an  identical  pair  of  equations  in  p'  and  kn.  Using  eqns  (29),  each  of  the  I-integrals 
can  be  reduced  to  a  one-dimensional  integral  over  s. 


To  complete  the  application  of  eqn  (27)  to  eqn  (22),  we  substitute  into  the 
remaining  integrals  involving  Ep  in  that  equation;  using  the  orthogonality  relations  in  eqns 
(18),  we  obtain  for  the  Nl"  order  approximate  admittance 


YNQ2Cn2(R/r) 
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+  jut. 


-2  -  «  -  M-  - 

“o  'oo  +  2“o  2  “m  'om  +  2  2  Vn 

m-1  m— 1  n-1 


mn 


(30) 


SIMPLIFICATION  AND  APPLICATION 


Eqn  (30)  is  homogeneous  in  the  fi's;  it  also  contains  many  appearances  of  factors 
(kir)Z1(kir)  and  (kjRJZjfkjR).  From  eqns  (7),  (6)  and  (18)  we  have 

C^rJZjCk,!-)  -  (Vklr>]'l(klr>(Jl<kir>Vkir)  -  J0<k,r>Vk,r>] 


(Yo<kir>)'1<kir) 
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'[Vkfr>] 

d(kj  r ) 


l[J0(kir)) 


Y0(kir)  dfkjr) 


The  expression  in  braces  is  the  Wronskian  of  the  zeroth-order  Bessel  functions  evaluated 
at  (kjr),  and  is  given  in  [8]  from  [9]  as  2/(*(kjr));  hence 


(kjr)21(kir)  -  2/[»Y0(kir)] 

and  simi larly 

(k.R)Z^ (k.R;  -  2/[*Y0(k.R)] 

Let  us  now  write 


(31) 


ai(k.r)Z1(k.r)  -  -  a.QQ  . 


Substituting  this  and  eqns  (31)  into  eqn  (30)  and  dividing  by  qq2  gives 
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where  the  integrals  are  now  the  limits  as  z'  tends  to  zero  of 
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According  to  the  asymptotic  formulae  for  Bessel's  Jq  (eg.  [9]),  the  factors  jQ(sr) 
and  Jq(sR)  tend  to  zero  like  s"i  as  s  tends  to  infinity;  hence  the  integrands  of  all  thse 
integrals  tend  to  zero  like  (exp(-z's)(s-^))  when  approaching  their  upper  limits.  This 
form  is  (algebraically)  convergent  even  when  z'  is  set  equal  to  zero,  so  we  can  now  take 
the  limit  and  obtain 
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The  integrands  of  these  integrals  have  several  false  singularities  and  (in  certain 
circumstances)  one  true  one.  We  shall  discuss  these  later;  first  we  complete  the  treatment 
of  eqn  (32),  which  is  now  an  easy  matter.  Taking  the  partial  derivatives  of  Yjq  with 
respect  to  each  of  the  o's  gives  us,  on  equating  them  to  zero  and  remembering  that 
Wi  =  Wi> 


a  r 

n  mn 


Vkmr> 


.W1 


0m 


(34) 


(m  -  1 . N) 


These  N  simultaneous  linear  equations  in  the  N  unknowns  am  determine  the  am  of  the 
N1*1  order  approximation,  and  the  corresponding  admittance  Yjq  is  then  obtained  from  eqn 
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(32).  One  further  simplification  is  possible;  if  we  multiply  each  of  eqns  (34)  by  its 
corresponding  om  and  sum  over  all  m,  we  obtain 
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and  so  eqn  (32)  reduces  to 
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(35) 


It  will  be  noticed  that  none  of  the  coefficients  or  right-hand  sides  in  these  equations 
depend  explicitly  on  N;  hence,  calculating  the  complete  set  for  a  given  N  provides  nearly 
all  the  information  needed  to  determine  the  approximate  admittances  for  every  smaller  N. 
This  is  very  useful,  as  it  allows  several  Y^  to  be  economically  calculated  and  then 
supplied  to  an  accelerated -convergence  routine  for  the  extrapolation  to  infinite  N. 

SINGULAR  BEHAVIOUR 

At  first  sight,  the  integrands  of  eqns  (33)  become  infinite  on  the  path  of 
integration,  at  s  =  0  in  I'oq,  s  =  km  in  lQm  and  at  s  =  km  or  kn  in  I'mn.  In  fact,  the 
numerators  of  the  integrands  have  compensating  zeros  at  all  these  points,  so  these 
singularities  are  illusory.  From  eqn  (5),  or  equivalently  from  eqns  (8),  (7)  and  (6), 
(Jo(sr)-(Yo(kmr)/Yo(kmR))Jo(sR))  has  a  zero  at  s  =  km  (and  since  Jq  is  an  even  function 
there  is  also  a  zero  at  s  =  -km,  although  this  is  not  directly  relevant);  similarly  there  is 
compensation  at  s  =  kn  in  I'mn,  even  when  kn  =  km.  However,  since  the  zeros  in  the 
numerators  of  the  integrands  are  generally  simple,  the  integrands  remain  finite  at  the  false 
singular  points  rather  than  vanishing  at  them;  this  behaviour  presents  a  numerical  problem 
which  is  discussed  further  below. 

There  is  no  such  problem  at  s  =  0.  Since  (Jo(sr)  -  Jq(sR))  is  of  order  (s2)  in 

that  neighbourhood,  the  integrand  of  J'qq  has  a  triple  zero,  rather  than  a  pole,  at  s  =  0; 

the  same  is  also  true  of  the  other  integrands,  although  again  this  is  not  directly  relevant 
since  they  do  not  have  apparent  poles. 

It  remains  to  consider  the  behaviour  of  the  integrands  near  s  =  kg.  Eqns  (33) 
were  established  with  kg  taken  as  strictly  complex,  while  s  is  purely  real  on  the  path  of 
integration;  so  in  principle  there  is  no  problem  here  either.  However,  when  the  medium 
in  region  B  is  lossless,  it  is  to  be  expected  that  the  behaviour  is  given  by  the  limit  of  the 
behaviour  with  a  slightly  lossy  medium  as  the  losses  tend  to  zero,  so  we  must  show  that 
this  limit  is  well-behaved  when  the  singular  point  at  s  =  kg  moves  onto  the  path  of 

integration. 

We  begin  by  confirming  that  the  singular  square  root  is  itself  well-defined.  The 

original  definition,  for  strictly  complex  kg,  was  that  (s2  -  kg2)i  tends  to  s  as  s  tends  to 
positive  real  infinity.  Over  the  real  half-range  of  s,  from  0  to  »,  (s2  -  kg2)  has  a  fixed 


positive  imaginary  part  and  a  real  part  which  is  large  and  positive  for  large  s,  falling  to 
zero  as  s  decreases  to  (real  part  of  (kg2))i  and  then  decreasing  further  to  -(real  part  of 
Accordingly,  the  phase  angle  of  (s2  -  kg2)i  starts  infinitesimally  greater  than 
zero  for  large  positive  s,  and  with  decreasing  s  it  increases  monotonically  to  (t/2  +  phase 
(kg)),  where  (-t/2)  <  phase  (kg)  <  0  since  kg  has  a  positive  real  part  and  a  negative 
imaginary  one.  In  the  limiting  case,  when  the  negative  imaginary  part  of  kg  tends  to 
zero,  the  phase  of  (s2  -  kg2)^  is  zero  for  s  >  kg  and  jumps  to  t/2  for  s  <  kg, 
corresponding  to  a  pure  positive  real  for  large  s  and  a  pure  positive  imaginary  for  small 
s. 


I 


The  easiest  way  to  prove  now  that  the  integrals  converge  for  real  kg  is  by  means 
of  the  transformation  previously  used  which  removes  the  singularity  for  all  kg.  Recalling 
the  treatment  of  eqn  (24)  above,  let  us  write  C  =  (s2  -  kg2)i,  where  the  definition  of 
the  square  root  has  just  been  examined;  then,  for  instance,  eqn  (33a)  gives 


‘sV°>!  t^y^)  -  jo[R/g2+ 4 )]' 

j  B  +  kB> 

7  (Vsr)  -  VsR>]2 

+  I  - -  -  ds 

§  s(s  -  kg)* 


I ' 

'oo 


dC 


(36) 


where:  S  is  some  convenient  fixed  large  positive  real  (greater  than  (real  part  of  (kg2))$), 
the  square  roots  in  the  Bessel  function  arguments  are  trouble-free  because  Jq  is  an  even 
function,  and  the  path  of  the  first  integral  passes  through  Q  =  0  when  kg  is  real.  (As  in 
treating  eqn  (24),  we  use  the  differential  transformation  s  ds  =  CdC.)  The  second  integral 
here  has  already  been  shown  to  converge,  having  satisfactory  behaviour  as  s  tends  to 
infinity.  The  integrand  of  the  first  integral  is  well-behaved  over  the  entire  finite  complex 
Q  plane,  the  apparent  poles  at  Q  =  zjkg  being  cancelled  by  (double)  zeros  of  the 
numerator,  and  the  range  of  integration  is  finite;  hence  this  integral  is  also  finite.  The 
other  I-integrals  can  be  treated  in  the  same  way  with  the  same  result;  their  apparent 
difficulties  at  s  =  ±km  become  apparent  difficulties  at  Q  =  t(km2  -  kg2)i  and  are  still 
compensated  as  discussed  earlier  in  this  section. 

NUMERICAL  EVALUATION 

We  must  now  consider  how  the  I-integrals  are  to  be  evaluated  in  practice. 

It  is  tempting  to  try  to  reduce  them  using  contour  integration,  but  there  does  not 
seem  to  be  any  satisfactory  way  of  doing  this.  There  is  an  extensive  literature  on  the 
Sommerfeld  problem  already  mentioned  (see  [14])  which  leads  to  integrals  of  somewhat 
similar  form;  but  all  the  techniques  employed  with  them  succeed  essentially  because  their 
integrands  depend  linearly  on  Bessel's  Jq,  and  not  on  Jq2  as  in  the  present  problem. 
We  must  therefore  use  an  essentially  numerical  approach. 


We  consider  first  the  problem  of  the  infinite  upper  limit  of  integration  in  eqns  (33). 
The  parts  of  the  integrands  involving  Bessel  functions  may  ail  be  written  as 


(vsr>  -  Vo<sR))(VsR)  -  Vo“R>) 
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where:  m  and  n  are  the  indices  used  in  eqns  (33),  yg  is  1,  and  ym  =  Yo(kmr)/Yg(kmR) 
for  non-zero  m.  For  large  s,  this  expression  may  be  reduced  using  the  asymptotic 
expansion  of  Jq  from  [11];  thus 
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The  remaining  factors  in  the  integrands  of  eqns  (33)  may  be  written  using  the  first-order 
binomial  theorem  as 
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where  kg  =  0.  At  first  sight  this  expression  is  inconsistent  with  eqn  (37),  in  which  terms 
of  two  inverse  orders  above  the  principal  terms  were  neglected.  However,  the 
development  of  eqn  (37)  proceeds  by  powers  of  (sr)-^  and  (sR)-1.  These  are  of 

comparable  magnitude,  since  R  is  roughly  (2r)  for  most  practical  cases,  and  (from  [8])  the 
first  mode  constant  kj  is  approximately  »/(R  -  r);  so  (sr)“*  (which  is  larger  than  (sR)-1 
and  therefore  determines  the  rate  of  decrease  of  successive  terms)  is  approximately 

k]/(rL),  where  L  is  the  lowest  value  of  s  to  which  the  asymptotic  form  is  applied.  We 
are  therefore  neglecting  in  eqn  (37)  terms  of  relative  size  (kj/(xL))2  but  retaining  terms 
of  relative  size  (($ks^  +  km*  +  kn2)/L^),  This  is  quite  reasonable;  at  the  N^1  order  of 
approximation  to  the  admittance,  km  and  kn  run  up  to  k]q,  which  is  roughly  (Nkj),  so 
the  ratio  of  the  two  relative  sizes  is  roughly  (*N)2  -  a  substantial  number  even  for 
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typical  N  of  6  or  8.  Also,  kg^  can  ^  much  larger  than  for  high-permittivity 
dielectrics,  and  moreover  the  imaginary  part  of  kg^  affects  the  real  part  of  the  admittance 
(which  is  usually  smaller  than  the  imaginary  part  and  more  difficult  to  calculate 
accurately).  It  is  therefore  quite  reasonable  to  use  expression  (39)  with  eqn  (38);  and  on 
combining  them  and  discarding  some  higher-order  product  terms,  we  obtain  as  the 
asymptotic  form  of  the  integrands  of  eqns  (33)  the  expression 
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The  integration  of  expression  (40)  between  the  lower  limit  L  (already  introduced)  and  the 
infinite  upper  limit  is  easily  performed,  using  asymptotic  integration  for  the  trigonometric 
terms;  eg. 
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where  some  higher  terms  have  been  neglected.  The  orders  of  approximation  in  eqn  (41) 
have  been  chosen  to  give  satisfactory  relative  accuracy  (a  few  parts  per  million)  for  N 
about  6  or  8,  R  comparable  with  (2r),  L  about  kgp  (note  that  L  is  not  dimensionless)  and 
the  magnitude  of  rg  (relative)  less  than  about  1000  at  the  highest  frequency  at  which  the 
coaxial  transmission -line  mode  propagates  alone  (the  cut-off  frequencies  of  the  various 
kinds  of  higher  mode  are  discussed  in  [8]).  (The  magnitude  of  the  smallest  of  the 
complete  integrals  in  eqns  (33)  is  readily  assessed  to  be  0(1 /(k^R)),  and  this  quantity  is 
used  as  the  reference.) 

To  complete  the  integrations  in  eqns  (33),  we  must  also  integrate  from  0  to  L. 
There  are  three  problems  associated  with  integration  over  this  range:  the  integrands 
oscillate  rapidly  (the  shortest  period  in  s  is  (*/R),  as  is  shown  by  the  asymptotic  analysis 
above,  and  this  is  of  the  order  of  kj  while  L  is  about  k$o)'.  m°st  of  the  integrands 
reduce  to  0/0  at  s  =  km  and/or  s  =  kn;  and  the  integrands  all  become  very  large,  or 
even  infinite,  when  s^  equals  the  real  part  of  kg^.  We  will  discuss  each  of  these 
problems  separately. 
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The  reduction  to  0/0  is  eliminated  as  follows.  Using  I'Hopital's  theorem,  we  can 
evaluate  each  integrand  accurately  at  the  false  singular  point  (or  points)  in  terms  of 
Bessel's  Jj.  Each  integrand  can  also  be  directly  evaluated  at  fixed  points  close  to  the 

false  singular  point,  with  a  modest  loss  of  significant  figures  (say,  3,  if  s  =  1.001  km  or 
0.999  km).  Within  the  range  spanned  by  the  fixed  points  we  can  then  use  interpolation 
which  keeps  the  loss  of  significant  figures  to  the  level  already  accepted  (say,  3)  provided 
the  order  of  interpolation  is  chosen  with  regard  to  the  accuracy  expected  in  the  Bessel 
functions.  In  the  author's  program  four  points  are  used  and  differences  up  to  the  third 
are  retained,  corresponding  to  12  significant  figures  in  the  raw  Bessel  functions. 

In  order  to  treat  the  other  two  problems  (the  rapid  oscillations  of  the  integrands 
and  the  behaviour  when  s2  is  near  kB2)  we  must  first  consider  how  the  numerical 

integration  should  be  done.  The  integrands  in  eqns  (33)  are  rather  complicated  functions 
of  the  integration  variable  s,  which  makes  it  desirable  to  keep  to  a  minimum  the  number 
of  evaluations  of  each  integrand;  accordingly,  Gaussian  integration  was  chosen,  using  a 
four-point  routine  for  eighth-order  accuracy.  An  elementary  integrator  of  this  kind  is 
normally  embedded  in  a  controlling  subroutine  which  calls  it  recursively  (or  in  a  loop), 
with  successively  smaller  sub-ranges  of  integration  until  the  desired  accuracy  is  achieved. 
In  the  present  problem,  we  can  eliminate  this  sub-division  by  making  a  virtue  out  of 
necessity;  since  the  shortest  period  in  s  is  (*/R),  it  is  convenient  and  satisfactory  to  use 
preset  sub-ranges  of  length  (ir/(4R)),  each  of  which  can  be  treated  to  satisfactory  accuracy 
using  a  single  four -point  elementary  Gaussian  integration.  This  method  has  the  happy 
consequence  that  all  the  integrands  in  eqns  (33)  -  ((N+l  )(N+2)/2)  of  them,  for  the  Nl” 
approximation  -  are  evaluated  at  the  same  values  of  s,  so  the  values  of  J o(sr)  and  Jq(sR) 

which  appear  in  them  can  be  determined  once  and  for  all;  this  leads  to  a  very  substantial 

saving  of  computation  time,  since  these  Bessel  functions  dominate  the  time  required  to 
evaluate  an  integrand. 

It  is  easy  to  identify  the  preset  sub-range  within  which  s2  becomes  equal  to  the 
real  part  of  kg2.  This  sub-range,  its  immediate  neighbour  above,  and  its  immediate 
neighbour  below  (if  the  frequency  is  high  enough  for  it  to  have  a  neighbour  below),  are 
collectively  treated  separately  from  all  the  other  sub-ranges,  thereby  isolating  the 
neighbourhood  of  the  singular  or  near-singular  point;  this  amount  of  isolation  is  sufficient 
to  keep  the  magnitudes  of  the  integrands  within  sensible  bounds  over  the  rest  of  the  range 
0  to  L.  Integration  over  the  three  special  sub-ranges  is  performed  using  the  elementary 
four-point  Gaussian  integrator  adaptively  (that  is,  with  recursively-nested  calls);  this  is 
slow,  but  is  accurate  and  straightforward  to  program.  The  three  sub-ranges  are  treated 
together  as  one,  but  the  first  of  the  internal  adaptive  subdivisions  is  always  taken  at  the 
point  s  =  (real  part  of  (kg2))i,  which  should  guarantee  that  no  attempt  will  be  made  to 
evaluate  an  infinite  integrand  when  kg  is  purely  real;  for  further  security  the  program 
changes  a  purely  real  kg  to  one  with  a  very  small  (negative)  imaginary  part. 

A  Gaussian  integrator,  whether  adaptive  or  elementary,  is  normally  designed  to 
integrate  a  real  function  over  a  real  range.  It  is  a  trivial  matter,  however,  to  extend  this 
to  the  integration  of  a  complex  function  of  a  real  variable  over  a  real  range,  since  all  the 
calculations  performed  during  the  integration  are  linear  in  the  calculated  values  of  the 
function.  The  integrands  in  eqns  (33)  are  only  complex  on  account  of  the  factor  (s2  - 
kg2)i  -  which  is,  of  course,  complex  for  some  values  of  s  even  if  kg  is  purely  real  - 
and  this  complex  square  root  can  be  evaluated  and  stored,  like  the  values  of  Jo(sr)  and 
Jq(sR),  for  each  of  the  repeated  values  of  s  as  described  above.  (A  complex  square  root 
can  be  extracted  using  two  real  square  roots.) 

The  rest  of  the  numerical  mathematics  is  comparatively  trivial  (see  (8)).  For  some 
convenient  moderately-sized  N  (6  or  8,  say)  the  I-integrals  are  calculated  and  their  values 
used  to  determine  the  os's  in  eqn  (34),  from  which  Yjq  follows  using  eqn  (35).  Using 
only  eqns  (34)  and  (35),  without  calculating  any  more  I-integrals,  we  can  also  determine 
Yj  for  i  =  1  (or  even  0)  to  i  =  (N-l).  The  resulting  set  of  Y's  can  then  be  processed 
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by  an  accelerated -convergence  routine  (such  as  the  one  described  in  [8])  to  find  the  true 
coaxial  terminating  admittance  Y. 

CONCLUSION 

A  rigorous  theoretical  treatment  of  the  open-ended  coaxial  line  with  infinite  flange 
has  been  presented,  together  with  all  the  essential  computational  details  required  to 
calculate  its  equivalent  admittance  at  the  plane  of  the  open  end.  The  treatment  requires 
the  numerical  evaluation  of  one-dimensional  integrals  only;  there  are  no  multiple  integrals 
involved. 


! 
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